---
title: "Replicate Main and SM analysis for US Global Health Aid Policy and Family Planning in Sub-Saharan Africa"
output: html_document
editor_options: 
  chunk_output_type: console
---

This RMarkdown file runs all of the analysis to reproduce the main and supplementary tables and figures for "US Global Health Aid Policy and Family Planning in Sub-Saharan Africa." Note that you must run the 3 data cleaning scripts first, ensure your directory structure follows the guidelines in the README, have all the packages in `setup_plgha()` installed, and run the `setup_plgha()` for this code to execute smoothly.

# Set-up
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

# load libraries and source functions 
source(here::here("scripts/functions.R")) 
setup_plgha(tidylog = F)

```


# Main Tables and Figures

## Tables 1 & 2: Summary of outcomes by exposure and timing of PLGHA

```{r main-tables}
sdp <- here("data-clean/plgha_sdp_df.rds") %>% 
  readRDS() %>% 
  filter(year >= 2014 & year <= 2019) %>%
  # Selections within packs
  mutate(
    services_prov = services_prov %>% 
      select(postabserprov, ancprov, pncprov, offerhivserv),
    mobile = mobile %>% select(recoutreach), 
    chv = chv %>% select(core8, long, short),
    across(
      c(fpprov, fpcharge, fprecout, fpoutday),
      ~.x %>% select(matches("core|long|short|emrg")) # drops "short" 
    )
  ) %>% 
  # What to include
  select(exposure, plgha_on, fpprov, fpcharge, fprecout, fpoutday, 
         services_prov, chv, mobile, facilityid) %>%
  # Unpack 
  unpack(where(is_tibble), names_sep = "_") %>%
  select(-c(fpprov_shortc, fpprov_longc, fpprov_shortnumc, fpprov_longnumc,
            fpprov_emrgc, fpprov_core8numc))

wmn <- here("data-clean/plgha_woman_df.Rds") %>% 
    readRDS() %>%
    # select only analytic years
    filter(year >= 2014 & year <= 2019)

dat <- list(
  wmn = wmn %>%
    select(exposure, plgha_on, fpcurruse, fptrad, fpemerg, fpstart, fpstartlarc,
           fpstop, pregnant, birth1yr, n_wmn),
  sdp = sdp %>% 
    select(-ends_with("num")) %>% 
    labelled::set_variable_labels(facilityid = "Number of facilities")
)

# Putting tables for `wmn` and `sdp` in a list from the jump makes it easier to repeat the same functions on both, then combine later.
dat_tabs <- dat %>% 
  map(
    ~.x %>% 
      mutate(
        stat_col = case_when( 
          exposure == "Low" & plgha_on ~ "low_on",
          exposure == "Low" & !plgha_on ~ "low_off",
          exposure == "High" & plgha_on ~ "high_on",
          exposure == "High" & !plgha_on ~ "high_off"
        )
      ) %>% 
      select(-c(exposure, plgha_on)) %>% 
      tbl_summary(
        by = stat_col,                          # creates four columns
        type = list(
          ends_with("num") ~ "continuous"
          ),
        statistic = list(
          all_continuous() ~ "{mean} ({sd})",   # mean and sd 
          all_categorical() ~ "{p}",            # percent 
          any_of("n_wmn") ~ "{n}",              # count cases 
          any_of("facilityid") ~ "{n_distinct}" # sdps appear in multiple rounds
        ),
        missing = "no"                          # drop counts of NA values
      ) %>% 
      add_stat_label() %>%  
      modify_table_body(
        ~.x %>%                                 # Hide stat_label for counts
          mutate(
            stat_label = ifelse(str_detect(stat_label, "^n"), NA, stat_label)
          )
      ) %>% 
      add_overall(last = TRUE) 
  )


# clean up formatting
dat_tabs <- dat_tabs %>% 
  map(function(tbl){
    tbl %>% 
      modify_table_body(~{
        out = .x %>%  
          mutate(
            var_label = var_label %>% str_remove(":.*"),
            label = label %>% str_remove(".*: ")
          ) %>% 
          mutate(
            row_type = case_when(
              var_label != label ~ "level",
              T ~ "label")
          ) 
        out %>% 
          group_by(var_label) %>% 
          slice(1) %>% 
          filter(row_type == "level") %>% 
          mutate(
            across(
              starts_with("stat") & !stat_label, 
              ~case_when(FALSE ~ .x)
            ),
            row_type = "label",
            label = var_label
          ) %>% 
          ungroup() %>% 
          bind_rows(out) %>% 
          mutate(variable = variable %>% str_remove("_.*")) %>% 
          arrange(var_label) %>% 
          mutate(stat_label = case_when(row_type == "label" ~ stat_label))
      }) %>% 
      gtsummary::italicize_levels()
  })

dat_tabs <- dat_tabs %>% 
  map(function(tbl){
    tbl %>% 
   # tbl_stack() %>% 
    modify_spanning_header(update = list(
      stat_1 ~ "PLGHA OFF (Obama Administration)\n2014-2016",
      stat_2 ~ "PLGHA ON (Trump Administration)\n2017-2019",
      stat_3 ~ "PLGHA OFF (Obama Administration)\n2014-2016",
      stat_4 ~ "PLGHA ON (Trump Administration)\n2017-2019"
    )) %>% 
    modify_header(update = list(
      label ~ "",
      stat_1 ~ "High Exposure",
      stat_2 ~ "High Exposure",
      stat_3 ~ "Low Exposure",
      stat_4 ~ "Low Exposure",
      stat_0 ~ "Overall"
    )) %>% 
    modify_footnote(update = everything() ~ NA) 
  })

# Arrange columns 
dat_tabs <- dat_tabs %>% 
   map(function(tbl){
     tbl %>% 
       modify_table_body(
        ~.x %>% 
          # PLGHA ON, then PLGHA OFF, then OVERALL 
         # relocate(stat_1, .before = last_col()) %>%
          relocate(stat_3, .after = stat_1) %>% 
          relocate(stat_0, .after = last_col())
      )
   })

# Arrange and drop unwanted rows     
dat_tabs <- dat_tabs %>% 
  map(function(tbl){
     tbl %>% 
        modify_table_body(
          ~.x %>% 
            # this mutate doesn't work anymore
            mutate(variable = variable %>% fct_relevel(
             "fpcurruse", "fptrad", "fpemerg", "fpstart", "fpstartlarc",
             "fpstop", "pregnant", "birth1yr", "fpprov", "fpcharge", "fprecout",
             "fpoutday", "services","mobile", "chv",   "n", "facilityid"
            )) %>% 
            arrange(variable)
        )
  })


# Export tables

## Table 1: Summary of facility-level outcomes by exposure and timing of PLGHA
sdp_xtable <- dat_tabs$sdp %>% 
  as_tibble() %>% 
  clean_names() %>%
  mutate(
    x = str_remove_all(x, "_"),
    x = str_replace(x, "%", "\\\\%"),
    x = str_replace(x, "&", "\\\\&"),
    x = str_replace(x, "core", "common"),
    x = case_when(
      !(str_detect(x, "\\\\%")) & !(str_detect(x, "Days out")) ~ paste0("\\hspace{3mm} \\textit{", x , "}"),
      T ~ x
    )
  ) %>%
  xtable(align = "llccccc", digits = 0)

autoformat(sdp_xtable)
addtorow <- list()
addtorow$pos <- list(0,0,0,30)
addtorow$command <- c(
  "& \\multicolumn{2}{c}{PLGHA OFF (2014-2016)} & \\multicolumn{2}{c}{PLGHA ON (2017-2019)} \\\\\n",
  " \\cmidrule(lr){2-3} \\cmidrule(lr){4-5} \n",
  "& High Exposure & Low Exposure & High Exposure & Low Exposure & Overall \\\\\n",
  "\\midrule \n"
)

print(
  sdp_xtable, 
  file = here("output/table_1.tex"),
  include.colnames = F, include.rownames = F, 
  booktabs = T, floating = F,
  add.to.row = addtorow, 
  sanitize.text.function = function(x) {x},
  format.args = list(big.mark = ","))


## Table 2: Summary of woman-level outcomes by exposure and timing of PLGHA
wmn_xtable <- dat_tabs$wmn %>% 
  as_tibble() %>% 
  xtable(align = "lLccccc", digits = 0)

autoformat(wmn_xtable)
addtorow <- list()
addtorow$pos <- list(0,0,0,8)
addtorow$command <- c(
  "& \\multicolumn{2}{c}{PLGHA OFF (2014-2016)} & \\multicolumn{2}{c}{PLGHA ON (2017-2019)} \\\\\n",
  " \\cmidrule(lr){2-3} \\cmidrule(lr){4-5} \n",
  "& High Exposure & Low Exposure & High Exposure & Low Exposure & Overall \\\\\n",
  "\\midrule \n"
)

print(
  wmn_xtable, 
  file = here("output/table_2.tex"), 
  include.colnames = F, include.rownames = F, 
  booktabs = T, floating = F,
  add.to.row = addtorow, tabular.environment = "tabularx",
  width="\\textwidth",
  format.args = list(big.mark = ","))


```

## Figure 1: Effect of PLGHA on Service Delivery

```{r run-sdp-models}
sdp <- here("data-clean/sdp.rds") %>% 
  readRDS() %>% 
  filter(year >= 2014 & year <= 2019)

# main models
sdp_did_results <-  sdp %>%
  mutate(services_prov = services_prov %>% select(
    postabserprov, ancprov, pncprov, offerhivserv
  )) %>%
  pack(fpoffered  = fpoffered) %>%
  sdp_did(
    clustervar = eaid, 
    fpprov, 
    fpcharge, 
    fprecout,
    fpoutday,
    services_prov,
    chv,
    mobile) 



# Select models and fix all the labels 
sdp_did_plots <- sdp_did_results %>%
  mutate(
    outcome = outcome %>% 
      str_replace("fpany_nummo", "fpanynummo") %>% 
      str_replace("fpnew_nummo", "fpnewnummo") %>% 
      str_replace("services_prov", "servicesprov")
  ) %>% 
  separate(outcome, sep = "_", into = c("outcome", "method")) %>% 
  mutate(
    outcome = if_else(outcome == "mobile", "chv", outcome),
    label = label %>%
      str_remove(".*: ") %>% 
      str_remove(" \\(if provided.*") %>%
      str_replace("Any core methods", "Any method"),
    label = if_else(
      outcome == "chv" & !str_detect(label, "Mobile"),
      paste("CHVs provide", label %>% str_replace("Any", "any")),
      label
    ),
    label = label %>% 
      str_wrap(24) %>%
      str_replace_all("\n", "<br>") %>%
      str_replace(
        "Post-abortion care",
        "Post-abortion care <span style='color:white;'>&&&&</span>"
      ),
    title = outcome %>% case_match(
      "fpprov" ~ "Family Planning Provision",
      "fpcharge"  ~ "Charges for Family Planning",
      "fprecout"  ~ "Recent Stockout of Family Planning",
      "fpoutday"  ~ "Length of Stockout in Days",
      "chv"  ~ "Community Health Volunteers and Mobile Outreach",
      "servicesprov" ~ "Other Health Service Provision"
    ),
    subtitle = outcome %>% case_match(
      c("fpcharge", "fprecout") ~ "Conditional on providing a given method",
      "fpoutday"  ~ "Conditional on being out of stock"
    )
  ) %>% 
  filter(
    model == "Controlled" & periodicity == "PLGHA",
    outcome %in% c("fpprov","fpcharge", "fprecout", "fpoutday", "chv",
                   "servicesprov"),
    method %in% c("core8", "long", "short", "emrg", "recoutreach") |
      outcome == "servicesprov"
  ) %>% 
  nest(.by = outcome)

# Set the plot order 
sdp_did_plots <- sdp_did_plots %>% 
  mutate(
    order = outcome %>% case_match(
      "fpprov" ~ "A",
      "fpcharge" ~ "B",
      "fprecout" ~ "C", 
      "fpoutday" ~ "D",
      "chv" ~ "E",
      "servicesprov" ~ "F"
    )
  ) %>% 
  arrange(order)

# Make a list of plots and arrange them in an array 
sdp_did_plots$plots <- sdp_did_plots$data %>%
  set_names(sdp_did_plots$outcome) %>%
  map(~{
    xmax <- .x %>% pivot_longer(starts_with("conf")) %>%
      pull(value) %>%
      abs() %>%
      max()
    xmin <- -xmax
    title <- unique(.x$title)
    subtitle <- unique(.x$subtitle)
    .x %>%
      did_plot() +
      labs(title = title, subtitle = if_else(!is.na(subtitle), subtitle, "")) +
      scale_color_manual(values = "black") +
      scale_y_continuous(limits = c(xmin, xmax), n.breaks = 5) +
      geom_text(aes(
        x = label,
        y = estimate,
        color = NULL,
        group = NULL,
        label = paste0(round(estimate, 3), " (", round(std.error, 3), ")")),
        size = 2.75,
        position = position_nudge(x = 0.2)) +
      theme(
        legend.position = "none",
        plot.title.position = "plot",
        plot.title = element_text(margin = margin(t = 18), size = 10),
        plot.subtitle = element_text(size = 9),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.x = element_blank(),
        axis.text.y = element_markdown(size = 8, hjust = 0),
        axis.text.x = element_text(size = 7),
        axis.title.x = element_text(size = 8)
      )
  }) 


cowplot::plot_grid(
  labels = "AUTO", 
  label_size = 9,
  plotlist = sdp_did_plots$plots, 
  ncol = 2
)

cowplot::ggsave2(
  filename = here("output/figure_1.pdf"), 
  height = 9, width = 7)

```

## Figure 2: Effect of PLGHA on Women

```{r run-wmn-models}
wmn <- here("data-clean/plgha_woman_df.Rds") %>% 
  readRDS() %>%
  filter(year >= 2014 & year <= 2019)


wmn_did_results <- wmn %>% 
    wmn_did(
      clustervar = eaid,
      fpcurruse, fptrad, fpemerg, fpstart, fpstop, pregnant, fpstartlarc, birth1yr,
      fpmodern, fplarc, fppill, fpinj, fpimp, fpcon, fpfc, fpiud
    ) 


fig_2 <- wmn_did_results %>%
  filter(model == "Controlled" & periodicity == "PLGHA") %>%
  filter(outcome %notin% c("fplarc", "fpinj", "fpmodern", "fppill",
                           "fpimp", "fpcon", "fpfc", "fpiud")) %>%
  mutate(
    label = as_factor(label),
    order = case_when(
      outcome == "fpcurruse" ~ 1, 
      outcome == "fptrad" ~ 2, 
      outcome == "fpemerg" ~ 3, 
      outcome == "fpstart" ~ 4, 
      outcome == "fpstartlarc" ~ 5, 
      outcome == "fpstop" ~ 6, 
      outcome == "pregnant" ~ 7, 
      T ~ 8
    ),
    label = fct_reorder(label, order)
  ) %>%
  did_plot() +
  geom_text(aes(x = label, y = estimate, 
                color = model, group = model,
                label = paste0(round(estimate, 3), " (", round(std.error, 3), ")")),
            position = position_nudge(x = 0.3)) +
  scale_color_manual(values = "black") +
  theme(legend.position = "none",
        panel.grid.major.y = element_blank(),
        axis.text.y = element_text(size = 12, hjust = 0)) 

ggsave(
  plot = fig_2,
  filename = here("output/figure_2.pdf"), 
  height = 4, width = 7)

```

# Supplementary Figure and Tables
## Fig S2: Disaggregation of Facility Results by Type of Modern Method Provided
```{r figure-s2}

fig_S2 <- sdp_did_results %>%
  filter(model == "Controlled" & periodicity == "PLGHA") %>%
  filter(outcome %in% c("fpprov_short", "fpprov_long", "fpprov_pill", 
                        "fpprov_injany", "fpprov_imp", "fpprov_iud",
                        "fpprov_con", "fpprov_fc")) %>%
  mutate(
    label = str_remove(label, "Provides: "),
    label = case_when(
      str_detect(label, "core") ~ "Any method",
      T ~ label
    ),
    label = fct_rev(as_factor(label)),
  ) %>%
  did_plot() +
  geom_text(aes(
    x = label, y = estimate, 
    color = model, group = model,
    label = paste0(round(estimate, 3), " (", round(std.error, 3), ")")),
    position = position_nudge(x = 0.2)) +
  scale_color_manual(values = "black") + 
  scale_y_continuous(breaks = pretty_breaks(4)) +
  theme(legend.position = "none",
        plot.title.position = "plot",
        panel.grid.major.y = element_blank(),
        axis.text.y = element_text(size = 12, hjust = 0),
        axis.title.x = element_text(size = 11)) 

ggsave(
  plot = fig_S2,
  filename = here("output/fig_S2.pdf"), 
       width = 6, height = 4)

```


## Fig S3: Robustness of Facility Results to Alternative Specifications and Exposure Definitions
```{r figure-s3}

select_sdp_outcomes <-  c("fpprov_core8", "fpprov_long", "fpprov_emrg",
                          "fpcharge_core8", "fpcurout_core8", "fpcurout_long",
                          "chv_core8", "services_prov_postabserprov", 
                          "services_prov_ancprov"
)

## plot alternative  exposures
fig_S3 <- sdp_did_results %>%
  filter(model %in% c("Controlled", "Uncontrolled", "FPRH Aid",
                      "US Aid % of Domestic Health Spending", "Log Aid") & 
           periodicity == "PLGHA") %>%
  filter(outcome %in% select_sdp_outcomes) %>%
  group_by(outcome) %>%
  fill(label) %>%
  ungroup() %>%
  mutate(
    model = as_factor(model),
    model = fct_relevel(model, "Controlled", "Uncontrolled", "FPRH Aid",
                        "US Aid % of Domestic Health Spending"),
    model = fct_recode(model,"Controlled (Base Model)" =  "Controlled"),
    model = fct_recode(model,"Ln(per capita GH aid)" =  "Log Aid"),
    label = as_factor(label)
  ) %>%
  did_plot() +
  theme(
    axis.text.y = element_text(size = 12, hjust = 0)
  )

ggsave(
  plot = fig_S3,
  filename = here("output/fig_S3.pdf"), 
  height = 7, width = 9.5)

```


## Fig S4: Robustness of Facility Results to Alternative PLGHA Timing
```{r figure-s4}

sdp_alt_timing_results <- map_dfr(
  list("plgha_on", "plgha_on_elect", "plgha_on_jan17", "plgha_on_jan18"),
  ~sdp %>% 
    filter(year >= 2014 & year <= 2019) %>% 
    mutate(services_prov = services_prov %>% select(
      postabserprov, ancprov, pncprov, offerhivserv
    )) %>%
    pack(fpoffered  = fpoffered) %>%
    mutate(plgha_on = !!sym(.x)) %>% # updates plgha_on with alt version
    sdp_did(
      clustervar = eaid, 
      fpprov, 
      fpcharge, 
      fpcurout,
      fpoutday,
      services_prov,
      chv,
      mobile,
      basic = T
    )  %>% 
    filter(model == "Controlled") %>% 
    mutate(model = .x) # names model according to alt plgha_on
)


sdp_drop2017_results <- sdp %>% 
  filter(year >= 2014 & year <= 2019 & year !=2017) %>% 
  mutate(services_prov = services_prov %>% select(
    postabserprov, ancprov, pncprov, offerhivserv
  )) %>%
  pack(fpoffered  = fpoffered) %>%
  sdp_did(
    clustervar = eaid, 
    fpprov, 
    fpcharge, 
    fpcurout,
    fpoutday,
    services_prov,
    chv,
    #mobile,
    basic = T
  )  %>% 
  filter(model == "Controlled") %>% 
  filter(outcome %in% select_sdp_outcomes)  %>%
  mutate(model = "Dropping 2017")


sdp_alt_timing_results <- sdp_alt_timing_results %>%
  filter(outcome %in% select_sdp_outcomes) %>%
  mutate(
    model = case_when(
      model == "plgha_on" ~ "Base Model (May 2017)",
      model == "plgha_on_jan17" ~ "Jan 2017",
      model == "plgha_on_jan18" ~ "Jan 2018",
      model == "plgha_on_elect" ~ "Nov 2016",
      T ~ model
    )
  ) %>%
  full_join(sdp_drop2017_results) 

fig_S4 <- sdp_alt_timing_results %>%
  did_plot() +
  theme(
    axis.text.y = element_text(size = 12, hjust = 0)
  )

ggsave(
  plot = fig_S4,
  filename = here("output/fig_S4.pdf"), 
  height = 7, width = 9.5)

```



## Fig S5: Robustness of SDP Result on Provision of Any Method of Family Planning to Dropping Individual Countries
```{r figure-s5}

base_model_fpprov <- sdp_did_results %>%
  filter(model == "Controlled" & periodicity == "PLGHA" & outcome == "fpprov_core8") %>%
  mutate(label = "Base Result") %>%
  mutate(model = "Full Sample")

sdp_unpack <- sdp %>%  
  filter(year >=2014 & year <=2019) %>%
  unpack(where(is_tibble), names_sep = "_")

loo_fpprov <- sdp_loo(data = sdp_unpack, outcome = fpprov_core8) %>%
  mutate(model = "Excluding Country") %>%
  mutate(label = left_out) 

fig_S5 <- loo_fpprov %>%
  full_join(base_model_fpprov) %>%
  did_plot() +
  geom_hline(yintercept = 
               base_model_fpprov$estimate[base_model_fpprov$model == "Full Sample"],
             color = "grey40") +
  theme(legend.position = "none",
        axis.text.y = element_text(size = 12, hjust = 0))

ggsave(
  plot = fig_S5,
  filename = here("output/fig_S5.pdf"), 
  height = 6.5, width = 8.5)

```



## Fig S6: Disaggregation of Woman-level Results by Type of Modern Method Used
```{r figure-s6}

fig_S6 <- wmn_did_results %>%
    filter(model == "Controlled" & periodicity == "PLGHA") %>%
    filter(outcome %in% c("fpmodern", "fplarc", "fpiud", "fpimp", "fpinj", 
                          "fppill", "fpcon", "fpfc", "fpimp")) %>%
    mutate(
        label = as_factor(label),
        order = case_when(
            outcome == "fpmodern" ~ 1, 
            outcome == "fplarc" ~ 2,
            outcome == "fpiud" ~ 3,
            outcome == "fpimp" ~ 4,
            outcome == "fpinj" ~ 5,
            outcome == "fppill" ~ 6, 
            outcome == "fpcon" ~ 7,
            outcome == "fpfc" ~ 8
        ),
        label = fct_reorder(label, order)
    ) %>%
    did_plot() +
    geom_text(aes(x = label, y = estimate, 
                  color = model, group = model,
                  label = paste0(round(estimate, 3), " (", round(std.error, 3), ")")),
              position = position_nudge(x = 0.3)) +
    scale_color_manual(values = "black") +
    theme(legend.position = "none",
          panel.grid.minor.x = element_blank(),
          panel.grid.major.y = element_blank(),
          axis.text.y = element_text(size = 12, hjust = 0)) 

ggsave(
  plot = fig_S6,
  filename = here("output/fig_S6.pdf"),
  width = 6, height = 4)

```



## Fig S7: Robustness of Woman-level Results to Alternative Specifications and Exposure Definitions
```{r figure-s7}


fig_S7 <- wmn_did_results %>%
  filter(model %in% c("Controlled", "Uncontrolled", "FPRH Aid",
                      "US Aid % of Domestic Health Spending", "Log Aid") & 
           periodicity == "PLGHA") %>%
  filter(outcome %notin% c("fplarc", "fpinj", "fpmodern", "fppill",
                           "fpimp", "fpcon", "fpfc", "fpiud")) %>%
  group_by(outcome) %>%
  fill(label) %>%
  ungroup() %>%
  mutate(
    model = as_factor(model),
    model = fct_relevel(model, "Controlled", "Uncontrolled", "FPRH Aid",
                      "US Aid % of Domestic Health Spending", "Log Aid"),
    model = fct_recode(model,"Controlled (Base Model)" =  "Controlled"),
    model = fct_recode(model,"Ln(per capita GH aid)" =  "Log Aid"),
    label = as_factor(label)
  ) %>%
  did_plot() +
  theme(legend.position = "bottom",
        axis.text.y = element_text(size = 12, hjust = 0)) 

ggsave(
  plot = fig_S7,
  filename = here("output/fig_S7.pdf"), 
  height = 7, width = 9.5)

```



## Fig S8: Robustness of Woman-level Results to Alternative PLGHA Timing

```{r figure-s8}

wmn_alt_timing_results <- map_dfr(
  list("plgha_on", "plgha_on_elect", "plgha_on_jan17", "plgha_on_jan18"),
  ~wmn %>% 
    filter(year >= 2014 & year <= 2019) %>% 
    mutate(plgha_on = !!sym(.x)) %>% # updates plgha_on with alt version
    wmn_did(
      clustervar = eaid,
      fpcurruse, fptrad, fpemerg, fpstart, fpstartlarc, fpstop, pregnant, birth1yr,
      basic = T
    ) %>% 
    filter(model == "Controlled") %>% 
    mutate(model = .x) # names model according to alt plgha_on
)

wmn_drop2017_results <- wmn %>% 
  filter(year >= 2014 & year <= 2019 & year !=2017) %>% 
  wmn_did(
      clustervar = eaid,
      fpcurruse, fptrad, fpemerg, fpstart, fpstartlarc, fpstop, pregnant, birth1yr,
      basic = T
    ) %>% 
  filter(model == "Controlled") %>%
  mutate(model = "Dropping 2017")

fig_S8 <- wmn_alt_timing_results %>%
  mutate(
    model = case_when(
      model == "plgha_on" ~ "Base Model (May 2017)",
      T ~ model
    )
  ) %>%
  full_join(wmn_drop2017_results) %>%
  mutate(
    label = as_factor(label),
    model = case_when(
      model == "plgha_on_jan17" ~ "Jan 2017",
      model == "plgha_on_jan18" ~ "Jan 2018", 
      model == "plgha_on_elect" ~ "Nov 2016",
      T ~ model
    )
  ) %>%
  did_plot() +
  labs(
    color = NULL
  ) +
  theme(legend.position = "bottom",
        axis.text.y = element_text(size = 12, hjust = 0)) 

ggsave(
  plot = fig_S8,
  filename = here("output/fig_S8.pdf"), 
  height = 6.5, width = 8.5)


```


## Fig S9: Robustness of Woman-level Result on Current Family Planning Use to Dropping Individual Countries
```{r figure-s9}

base_model_fpcurruse <- wmn_did_results %>%
  filter(model == "Controlled" & periodicity == "PLGHA" & outcome == "fpcurruse") %>%
  mutate(label = "Base Result") %>%
  mutate(model = "Full Sample")

loo_fpcurruse <- wmn_loo(data = wmn, outcome = fpcurruse) %>%
  mutate(model = "Excluding Country") %>%
  mutate(label = left_out) 

fig_S9 <- loo_fpcurruse %>%
  full_join(base_model_fpcurruse) %>%
  did_plot() +
    geom_hline(yintercept = 
                 base_model_fpcurruse$estimate[base_model_fpcurruse$model == "Full Sample"],
               color = "grey40") +
    theme(legend.position = "none", 
          axis.text.y = element_text(size = 12, hjust = 0)) 

ggsave(
  plot = fig_S9,
  filename = here("output/fig_S9.pdf"), 
  height = 6.5, width = 8.5)

```


## Fig S10: Sample Countries and Exposure to PLGHA
```{r figure-s10}

africa <- st_read("data-raw/Other/afr_adm.shp") %>%
  filter(!is.na(NAME)) %>%
  mutate(
    country = str_to_title(NAME, locale = "en"),
    country = case_when(
      country == "Zaire" ~ "Democratic Republic of the Congo",
      T ~ country
    )
  ) %>%
  st_set_crs(4326) 

exposure <- here("data-clean/exposure.rds") %>% 
  readRDS()


africa <- left_join(africa, exposure) 

fig_S10 <- africa %>%
  ggplot() +
    geom_sf(aes(fill = exposure)) +
    scale_fill_npg(
      na.translate = F,
      labels = c("Low (< median\nper capita aid)",
                 "High (> median\nper capita aid)") ) +
    labs(fill = "Exposure to PLGHA") +
    theme_bw() +
    theme(
      panel.grid = element_line(colour = "transparent"),
      panel.border = element_blank(),
      legend.text = element_text(size = 6),
      legend.title = element_text(size = 8),
      panel.background = element_blank(),
      plot.background = element_blank(),
      legend.background = element_blank(),
      axis.text = element_blank(),
      axis.ticks = element_blank(),
      legend.position = c(0.28, 0.3))

ggsave(
  plot = fig_S10,
  filename = here("output/fig_S10.pdf"), 
  width = 4, height = 4)


```


## Fig S11: Countries and PMA Survey Years, By Exposure
```{r figure-s11}

exp_year_tab <- wmn %>%
    count(country, exposure, year, name = "Women", plgha_on) 

update_geom_defaults("text", aes(size = 3))

fig_S11 <- exp_year_tab %>% 
    filter(year >= 2014 & year <= 2019) %>%
    count(country, exposure, year) %>% 
    mutate(
        country = country %>% fct_rev,
        exposure = exposure %>% paste0("\nExposure") %>% fct_rev) %>% 
    ggplot(aes(x = year, y = country, fill = exposure)) + 
    geom_tile(alpha = 1, height = 0.5, width = 0.8) + 
    geom_text(aes(label = year)) +
    facet_grid(rows = vars(exposure), scales = "free") + 
    theme_minimal() +
    theme(
      text = element_text(size = 12),
      axis.text.x = element_blank(),
      strip.text.y = element_text(angle = 0),
      panel.grid = element_blank(),
      legend.position = "none"
    ) + 
    scale_fill_npg() +
    scale_y_discrete(labels = ~str_wrap(.x, width = 20)) + 
    labs(x = NULL, y = NULL)


ggsave(
  plot = fig_S11,
  filename = here("output/fig_S11.pdf"), 
  width = 6, height = 4)

```

## Fig S12: Generalized difference-in-differences

```{r figure-s12}
sdp_gdid <- sdp_did_results %>%
  filter(model == "Controlled" & periodicity == "GDiD" & outcome == "fpprov_core8") %>%
  group_modify(~add_row(.x, .before = 0)) %>%
  mutate(
    period = case_when(
      is.na(period) ~ "2016",
      T ~ period),
    estimate = case_when(
      is.na(estimate) ~ 0,
      T ~ estimate)
  ) %>%
  fill(label, .direction = "up") %>%
  fill(model, .direction = "up") %>%
  ungroup() %>%
  mutate(
    label = str_replace(label, "Any core methods", "any FP method"),
    label = str_replace(label, "Provides:", "SDP provides"),
    label = str_remove(label, " \\s*\\([^\\)]+\\)"), # remove (if provided)
    label = as_factor(label),
    period = fct_rev(factor(period,
                               levels = c("2014", "2015",
                                          "2016", "2017",
                                          "2018", "2019")))
    ) %>%
  gdid_plot() +
  scale_color_manual(values = "black") +
  theme(legend.position = "none",
       plot.title.position = "plot",
       axis.text.y = element_text(size = 10)) 

wmn_gdid <- wmn_did_results %>%
  filter(model == "Controlled" & periodicity == "GDiD" & outcome == "fpcurruse") %>%
  group_modify(~add_row(.x, .before = 0)) %>%
  mutate(
    period = case_when(
      is.na(period) ~ "2016",
      T ~ period),
    estimate = case_when(
      is.na(estimate) ~ 0,
      T ~ estimate)
  ) %>%
  fill(label, .direction = "up") %>%
  fill(model, .direction = "up") %>%
  ungroup() %>%
  mutate(
    label = as_factor(label),
    period = fct_rev(factor(period,
                               levels = c("2014", "2015",
                                          "2016", "2017",
                                          "2018", "2019")))
    ) %>%
  gdid_plot() + scale_color_manual(values = "black") +
  theme(legend.position = "none") 



fig_S12 <- sdp_gdid + wmn_gdid +
  plot_annotation(tag_levels = "A")

ggsave(
  plot = fig_S12,
  filename = here("output/fig_S12.pdf"), 
  height = 4, width = 8.8)

```


# Supplementary Tables
## Table S2: Comparison of Exposure by Different Classifications
```{r table-s2}

aid <- here("data-clean/exposure.rds") %>% 
  read_rds()

#median_aid <- unique(aid$med_aid)

table_S2 <- aid %>% 
  dplyr::select(country, aidpc, exposure, us_pct_ghes, exposure_ghes,
                fprhpc, exposure_fprh) %>%
  arrange(-aidpc) %>%
  as_tibble() %>% 
  xtable(align = "llllllll", digits = c(0,0,2,0,1,0,3,0))


addtorow <- list()
addtorow$pos <- list(0,0,0)
addtorow$command <- c(
  "Country & \\multicolumn{3}{c}{US Bilateral Global Health Aid} & \\multicolumn{2}{c}{US \\% of Gov. Health Exp.}   & \\multicolumn{2}{c}{US bilateral FPRH Aid} \\\\\n",
  " \\cmidrule(lr){2-4} \\cmidrule(lr){5-6} \\cmidrule(lr){7-8}\n",
  "& Per Capita & Exposure & Per Capita & Exposure & Per Capita & Exposure \\\\\n"
)



print(
  table_S2, 
  file = here("output/table_S2.tex"), 
  include.colnames = F, include.rownames = F, 
  booktabs = T, floating = F,
  add.to.row = addtorow, 
  format.args = list(big.mark = ",")
)


```


